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ABSTRACT 

Thermodynamic flowline and plume models for the ice shelf-ocean system simplify the ice and ocean 
dynamics sufficiently to allow extensive exploration of parameters affecting ice-sheet stability while including 
key physical processes. Comparison between geophysically and laboratory-based treatments of ice-ocean 
interface thermodynamics shows reasonable agreement between calculated melt rates, except where steep 
basal slopes and relatively high ocean temperatures are present. Results are especially sensitive to the poorly 
known drag coefficient, highlighting the need for additional field experiments to constrain its value. These 
experiments also suggest that if the ice-ocean interface near the grounding line is steeper than some 
threshold, further steepening of the slope may drive higher entrainment that limits buoyancy, slowing the 
plume and reducing melting; if confirmed, this will provide a stabilizing feedback on ice sheets under some 
circumstances. 


1. Introduction 

Ice shelves, the floating extensions of the outlet gla- 
ciers and ice streams that drain inland ice sheets, cover 
approximately 40% of the Antarctic continental shelf 
(Williams et al. 1998). As the underlying ocean is ef- 
fectively isolated from any atmospheric influence, cir- 
culation beneath ice shelves is primarily thermohaline in 
nature, driven by melting and freezing at the shelf base, 
with tides also contributing to vertical mixing (MacAyeal 
1984). This local circulation is of global interest for two 
reasons. First, the outflow of cooled, freshened Ice Shelf 
Water (ISW) formed by sub-shelf melting contributes to 
the formation of Antarctic Bottom Water (AABW), a 
key driver of thermohaline circulation. Second, the pres- 
ence of warm, dense Circumpolar Deep Water (CDW) 
in the Amundsen Sea has led to high melt rates (tens of 
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meters per year), thinning and retreat of ice shelves in 
this region of West Antarctica, and acceleration of in- 
land ice streams (Rignot 1998; Rignot and Jacobs 2002; 
Payne et al. 2004; Shepherd et al. 2004). Although 
melting of floating ice has little direct effect on sea level 
(Jenkins and Holland 2007), loss of the buttressing pro- 
vided by ice shelves causes increased flux of grounded 
ice into the ocean, drawing down the interior ice sheet 
and contributing to sea level rise (e.g., Dupont and Alley 
2005, 2006). Uncertainty about this process was a key 
factor in the most recent Intergovernmental Panel on 
Climate Change (IPCC) report providing projections of 
sea level rise “excluding future rapid dynamical changes 
in ice flow” (Solomon et al. 2007). While recent com- 
munity modeling efforts such as Sea-level Response to 
Ice Sheet Evolution (SeaRISE; http://websrv.cs.umt.edu/ 
isis/index.php/SeaRISE_Assessment; Bindschadler et al. 
2013) and Ice2Sea (http://www.ice2sea.eu/) have provided 
more realistic bounds on glacially driven sea level rise, 
several studies (Walker et al. 2008; Gagliardini et al. 
2010) have shown that both the magnitude and spatial 
distribution of basal melting strongly affect ice shelf 
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buttressing, suggesting that ice-sheet models will require 
accurate oceanic forcing to project sea level rise with 
greater precision. 

However, modeling of sub-ice shelf circulation re- 
mains subject to considerable uncertainty, largely be- 
cause of the inaccessibility of these regions and the 
resulting scarcity of observations. Until the recent de- 
velopment and deployment of autonomous submersibles 
(Nicholls et al. 2006; Jenkins et al. 2010a, 2012), obser- 
vations were limited to ship-based measurements in the 
open ocean near ice fronts (Jacobs et al. 1996; Smethie 
and Jacobs 2005; Jacobs et al. 2011) and a handful of 
studies in which instruments were lowered through 
holes drilled in ice shelves (e.g., Nicholls et al. 2009). 
Several key elements in the modeling of sub-ice shelf 
circulation are borrowed from analogous oceanic phe- 
nomena that have been far more extensively observed. 
Dynamical models of buoyant ISW plumes (e.g., Jenkins 
1991; Holland and Feltham 2006) are based upon treat- 
ments of dense overflows down slopes in laboratory ex- 
periments (Ellison and Turner 1959, 1973) and in the ocean 
(e.g., Killworth 1977; Bo Pederson 1980; Stigebrandt 
1987; Arneborg et al. 2007). Causing greater uncertainty, 
theories of thermodynamics at the ice shelf-ocean inter- 
face (e.g., Holland and Jenkins 1999) rely upon the as- 
sumption that observations of melting and vertical mixing 
below sea ice (McPhee 1990, 1992; McPhee et al. 1999; 
McPhee 2008) are applicable on larger scales. While 
these sea ice observations are both extensive and highly 
detailed, there are not yet enough sub-ice shelf mea- 
surements to confirm this hypothesis, and it is uncertain 
whether basal roughness and ocean currents are suffi- 
ciently similar. Recent work by Jenkins et al. (2010b) 
shows that melt rates observed at one site beneath the 
Ronne Ice Shelf are consistent with multiple parame- 
terizations of the turbulent boundary layer, highlighting 
the need for further observations under a wide range of 
conditions to narrow the uncertainties. 

In this study, we examine several uncertainties in in- 
terface thermodynamics and their likely impact on ice 
dynamics. We begin by assessing whether three plausi- 
ble versions of the thermodynamical equations produce 
melt rates sufficiently different to significantly affect the 
ice-ocean system. Our analysis of the results suggests a 
significant role for the drag coefficient, motivating a 
sensitivity study on this parameter. We then apply 
a simplified, fully coupled ocean-ice shelf-ice stream 
model to assess the range of grounding-line retreat re- 
sulting from uncertainty in the drag coefficient. Finally, 
because ice shelf geometry strongly affects basal melting 
(Jenkins 1991; Walker and Holland 2007; Little et al. 
2009), we repeat some of our fixed-shelf experiments for 
multiple ice shelf depth profiles, in some cases finding an 


interesting relationship between melt rates and basal 
slope. 

2. Ocean model 

a. Ocean plume dynamics 

This study uses a slightly modified version of the one- 
dimensional ocean plume model introduced by Jenkins 
(1991): 

yr(UD) = e + m, (1) 

ax 

^ D) ,- cD (e^ gsine -c d u^ ( 2) 

j x (TUD) = T w e+T b m + (T b -T)y T , and (3) 

Y x ( SUD ) = V + V* + ( s „ ~ s )y s > ( 4 ) 

where U, D, T, S, and p indicate plume velocity, thickness, 
temperature, salinity, and density, respectively; p 0 is the 
reference density; g is the acceleration due to gravity; 9 is 
the slope of the ice shelf base; C lt is the drag coefficient; 
and T w , S w , and p w are the temperature, salinity, and 
density of the ambient seawater outside the plume, re- 
spectively. We modify the original momentum equation 
by including a geostrophic factor e in (2) to parameterize 
the otherwise neglected Coriolis force (cf. Wright and 
Stocker 1991; Walker et al. 2009; Parizek and Walker 
2010). While this model cannot fully capture the effect of 
water column thickness on depth-integrated barotropic 
flow, studies using three-dimensional isopycnic-coordinate 
models with explicit mixed layers (e.g., Little et al. 2009) 
have shown that basal slope strongly affects melt rates, 
providing some justification for the plume approxi- 
mation. We also modify the calculation of the entrain- 
ment rate 

e = e Q esm(8)U (5) 

by setting e 0 = 1.2 C]j 2 as in Stigebrandt (1987), rather 
than using constant e 0 as in the original model. This 
parameterization allows entrainment to depend on the 
drag coefficient, as in Bo Pederson (1980) and Arneborg 
et al. (2007), while remaining simple enough for use with 
a reduced model. 

b. Laboratory-based interface thermodynamics 

The melt rate m, ice shelf basal temperature and sa- 
linity T b and S h , and turbulent exchange velocities y T 
and y s are calculated by a separate set of equations for 
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thermodynamics at the ice-ocean interface. We typi- 
cally use a three-equation formulation (Hellmer and 
Olbers 1989; Holland and Jenkins 1999), in which heat 

and salt balances are calculated across a 
assumed to be at the local freezing point: 

thin sub-layer 

T b =\ l S B + A 2 + A 3 p B , 

(6) 

C pM T surt ~ T b ) + C V ^P T ~ T b ) = ™ L / 

, and (7) 

7,5 (-* ~ $ g) — mS B , 

(8) 


where A, are coefficients in the linearization of the 
freezing point (Millero 1978), p B is the pressure at the 
ice shelf base, c pi and c pw are the specific heat capacities 
of ice and seawater, Lf is the latent heat of fusion, and 
r surf is the ice shelf surface temperature. The turbulent 
exchange velocities are calculated using the equations 
derived by Jenkins (1991) from the analyses of labora- 
tory studies reported by Kader and Yaglom (1972, 1977): 
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where v w , Pr, and Sc are the kinematic viscosity and 
molecular Prandtl and Schmidt numbers of seawater. In 
subsequent experiments, we will refer to the combina- 
tion of (6)-(8) with (9) as “laboratory” three-equation 
thermodynamics. 


c. Geophysically based interface thermodynamics 

Several studies based on field observations have led 
to alternative formulations of ice-ocean interface ther- 
modynamics. While these methods lead to simpler sets 
of equations, their primary motivation is to describe in- 
terface thermodynamics in terms of variables that can be 
measured in a geophysical setting. 

Jenkins et al. (2010b) found that melt rates observed 
at the base of Ronne Ice Shelf by phase-sensitive radar 
could be accurately calculated by a suitably tuned ver- 
sion of (9), but that their data could be matched equally 
well using a simpler parameterization in which the 
Stanton numbers (C ] J 2 Vr, C 1 / 2 !^) are constants. Not- 
ing that the complexity of expressions like (9) has not 
been shown necessary on geophysical scales, and that in 
practice this formula gives nearly constant values because 


of the dominance of molecular diffusion in the interface 
sublayer, they recommended the use of constant Stanton 
numbers with (6)-(8). 

McPhee (1992) and McPhee et al. (1999) recom- 
mended (at least under sea ice) that (8) could be drop- 
ped by assuming the interface salinity to be equal to the 
plume salinity and parameterizing the rate limiting 
process of salt diffusion through the molecular sublayer 
by a suitable choice of effective heat transfer velocity, 
leaving the two-equation formulation 

T b = AjS + A 2 + A 3 /?g and (10) 

c piM T surf - T B ) + c pv/ y TS (T - T B ) = rhL f , (11) 

which agrees well with measurements of heat fluxes 
beneath sea ice of widely varying roughness. Jenkins 
et al. (2010b) found that this formulation could also 
match their observations when used with constant Stanton 
number (Cf 2 T re). although they cautioned that it is likely 
less applicable across a broader range of oceanographic 
conditions than three-equation thermodynamics. 


3. Ocean experiments 

Our experiments using the plume model beneath a 
fixed ice shelf emulate the study of Holland et al. (2008) 
using a three-dimensional isopycnic-coordinate ocean 
model. We consider shelves of lengths 275 and 550 km, 
with depth ranging from 600 m at the grounding line to 
200 m at the ice front according to the formula 


z = 


(a 7 + x) 


lln ’ 


( 12 ) 


where a t are chosen to give the desired depths at 
grounding line and ice front, and the exponent n allows 
the ice shelf profile (Fig. 1) to vary from linear ( n = — 1) 
to highly concave (n = 4). The grounding-line depth has 
been reduced from 1000 in in the earlier study to avoid 
situations involving low melt rates (at lower ocean tem- 
peratures) and steep basal slopes, which may cause the 
plume to reach neutral buoyancy and separate from the 
ice shelf base, terminating the simulation. This shallower 
grounding line allows us to run experiments with deep 
ocean temperatures ranging from — 1.8° to 2.0°C at 0.2°C 
intervals. We apply the temperature-salinity profiles 
used by Holland et al. (2008), setting the vertical salinity 
gradient and deep ocean temperature constant for depths 
greater than 200 m. 

In presenting these experiments, we should emphasize 
that, unlike Jenkins et al. (2010b), we do not have suf- 
ficient observations to evaluate the accuracy of each 
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Fig. 1. Depth profiles of ice shelf bases with length 275 km 
generated by (12) for n = 4 (steepest), 3, 2, 1, and —1 (linear). 


thermodynamical scheme. Instead, we will assess how 
closely the geophysically based two- and three-equation 
constant Stanton number methods and the laboratory- 
based three-equation formulation agree. We will also 
examine which parameters are primarily responsible for 
differences in model results, with an eye toward guiding 
further observations. 

a. Comparing geophysically and laboratory-based 
thermodynamics 

In our first set of experiments, all parameters are as 
determined by Jenkins et al. (2010b) through tuning the 
various methods to match their observations. The lab- 
oratory three-equation method uses C,j = 0.0062, while 
the geophysical methods use a somewhat higher value 
[C d = 0.0097, with Cy 2 T s = 3.1 X 10~ 5 , Cf-V T = 0.0011, 
and ( C 1 J 2 Tts ) = 5.9 X 10~ 4 ]. For the flattest shelf, the 
550-km linear profile, maximum and mean melt rates 
from the geophysical methods are 9%-ll% higher than 
laboratory. For the steepest shelf, the 275-km nonlinear 
profile with n = 4, the two-equation constant Stanton 
number method produces the lowest melt rates, while 
the three-equation constant Stanton method closely 
matches the laboratory method (Fig. 2). 

In an effort to separate the effects of thermodynamical 
scheme from those of parameterization, we also consider 
experiments in which the laboratory three-equation 
method uses the same drag coefficient ( C,/ = 0.0097) as 
the constant Stanton number methods. These runs pro- 
duce better agreement between methods for the 550-km 
linear ice shelf by increasing laboratory melt rates to less 
than 3% higher than either geophysical method. How- 
ever, agreement between methods worsens for the 275-km 
nonlinear shelf (Fig. 2, bottom). 




Fig. 2. Max melt rates for (top) 550-km linear (n = —1) and 
(bottom) 275-km nonlinear (n = 4) ice shelf profiles, using C& = 
0.0097 for constant Stanton number runs and both C <t = 0.0062 and 
Cj = 0.0097 for laboratory three-equation thermodynamics. Two- 
equation runs use (10) and (11), three-equation runs use (6)— (8), 
and laboratory-based runs also use (9) for exchange velocities. 
Note differing y-axis scales. 

The differences between the preceding sets of runs 
can be explained by examining the three roles that the 
drag coefficient C,/ plays in the model. First, for the 
laboratory three-equation method, C 1 / 2 is a factor in 
calculating the turbulent exchange velocities using (9), 
so that increasing the drag coefficient produces more 
vertical mixing (and thus more melting) at a given plume 
velocity. (For the other methods, the drag coefficient is 
already included as part of the constant Stanton num- 
bers.) Second, a factor of C)j 2 also appears in our en- 
trainment parameterization (5), so that at a given plume 
velocity, increasing the drag coefficient will increase 
entrainment of warmer, saltier ambient water into the 
plume, producing more melting. Third, the drag co- 
efficient appears in the momentum equation (2), where 
a higher value will decrease plume velocity, working 
against vertical mixing (at the ice-ocean interface) and 
entrainment (at the plume-ambient water interface). It 
is thus not immediately clear whether increasing the drag 
coefficient will increase or decrease basal melting, but we 
will find later that the influence of drag coefficient on 
entrainment is the strongest of these effects. 

For the 550-km linear shelf, the lower drag coefficient 
( C,i = 0.0062) laboratory runs produce a slightly faster- 
flowing plume than the (three equation) constant Stan- 
ton number runs. Despite this advantage in velocity, the 
lower drag coefficient leads to a smaller thermal Stanton 
number [mean C 1 / 2 I> = 9.7 X 10 1 from (9), versus 
constant value of 1.1 X 10 3 ] and a smaller thermal 
exchange velocity y T . The lower C,/ also leads to a 
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slightly lower entrainment rate from (5), resulting in a 
smaller temperature contrast T — Tb between the plume 
and the interface sublayer. As solving (6)-(8) for melt 
rate gives, to leading order, m «■ y r (T — Tb), the con- 
stant Stanton number runs produce approximately 10% 
greater melt. When the laboratory runs are repeated 
with C,/ = 0.0097, the higher drag coefficient more than 
offsets a decrease in velocity, leading to increased values 
of thermal exchange velocity and entrainment rate and 
approximately 2% greater melt than for the constant 
Stanton number runs. 

For the 275-km nonlinear shelf, the extra velocity 
gained by the laboratory (C,i = 0.0062) runs compen- 
sates for the lower drag coefficient, leading to thermal 
exchange velocities and entrainment rates only slightly 
lower than those of the constant Stanton number runs. 
The two methods thus agree rather closely for this case. 
In contrast to the linear shelf case, increasing the drag 
coefficient for the laboratory runs worsens the agree- 
ment significantly. The plume slows slightly, but a high 
thermal Stanton number leads to an increased turbulent 
exchange velocity. With the entrainment rate for the 
laboratory runs actually slightly exceeding that of the 
constant Stanton number runs, the laboratory runs now 
produce greater melting than the constant Stanton number 
runs. 

Overall, these results indicate that three-equation 
thermodynamics with constant Stanton numbers and the 
laboratory-based formulation match rather well, if cal- 
culated melt rates are of greatest interest. Two-equation 
thermodynamics can also provide a reasonable match, 
except in cases involving both high ocean temperatures 
and high ice shelf basal slopes (cf. Fig. 2). Differences 
between the two three-equation methods can likely be 
further reduced by suitable choices of model parameters, 
although determination of the drag coefficient presents its 
own difficulties. 

It must be noted that Jenkins et al. (2010b) urge 
caution in the use of their relatively high values for the 
drag coefficient, pointing out that their data do not allow 
independent evaluation of both drag and turbulent trans- 
fer coefficients. While their values for Stanton numbers 
are well constrained, deriving a value for the drag co- 
efficient requires the assumption that the turbulent 
transfer coefficient found by McPhee (1992) for sea ice is 
valid beneath ice shelves. Uncertainties in the measure- 
ment of temperature and velocity also contribute to un- 
certainty in the drag coefficient. 

b. Assessing sensitivity to drag coefficient 

The difficulty of precisely determining the drag co- 
efficient even where observations exist motivates sen- 
sitivity studies on this parameter. While we have seen 




Fig. 3. Max melt rates for (top) 550-km linear (n = —1) and 
(bottom) 275-km nonlinear (n = 4) ice shelf profiles, using labo- 
ratory three-equation thermodynamics. Note differing y-axis scales. 


that the three-equation methods are comparable, we use 
laboratory thermodynamics because evaluation of (9) is 
more straightforward than determining the proper con- 
stant Stanton numbers for each value of Cy, the extra 
computational expense is negligible for our one-dimensional 
model. We repeat the experiments already described 
with three commonly used values of Cj (0.0015, 0.0025, 
and 0.0035) and the two values suggested by Jenkins 
et al. (2010b) (0.0062 and 0.0097). Both the 550-km 
linear and 275-km nonlinear (Fig. 3) shelf profiles show 
significantly higher melt rates as C',/ increases, indicating 
that the effects of enhanced vertical mixing and en- 
trainment are stronger than the effect of reduced veloc- 
ity. Compared to C t / = 0.0025, which has been a standard 
estimate since MacAyeal (1985), Ca = 0.0097 produces 
36%-40% higher mean and maximum melt rates for the 
550-km linear shelf; for the 275-km nonlinear shelf, with 
faster plume flow, mean and maximum melt rates are up 
to 46% higher. Typical values of the diffusive Stanton 
number nearly double (from 2.20 X 10~ 5 to 4.33 X 10~ 5 ) 
and the entrainment rate increases by nearly half (from 
1.32 X 10- 5 to 1.93 X 10 _5 ms _1 for the 275-km non- 
linear shelf at 0°C); together, these effects more than 
make up for the plume slowing (from 8.34 to 6.05 cm s 1 
for the 275-km nonlinear shelf at 0°C). 

To test the relative importance of drag-dependent 
processes, we repeat the above experiments with the 
entrainment parameter e 0 set to the constant value used 
in Jenkins (1991). When entrainment no longer depends 
on the drag coefficient, we find that the effect of reduced 
velocity overcomes increased Stanton numbers and 
mean melt rates drop as C, t increases. The relative change 
is largest for linear ice shelves, which lack the steep near- 
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grounding-line slopes needed to establish plume velocity 
against higher drag coefficients. For the 275-km nonlinear 
ice shelf, maximum melt rates for Cd = 0.0097 can be up 
to 15% higher than for Cd = 0.0025, even though mean 
melt rates are lower. This effect is driven by steep near- 
grounding-line slopes that allow plume velocities fast 
enough to take advantage of the increase in Stanton 
numbers at high Cd', the relative difference is largest at 
low temperatures, where the slope is most important in 
driving melting, and decreases as the ocean warms. Still, 
the overall result is decreasing melt as the drag co- 
efficient is increased, the opposite of what is found in 
experiments with drag-dependent entrainment. We thus 
conclude that entrainment is the most significant of the 
drag-dependent processes in our model. 


T, 


T b u , h > Ilf 

0, h < h f ’ 


(16) 


where r y and T h are constant coefficients, and hf is the 
hydrostatically determined flotation thickness. At the 
downstream end (x = L x ), we impose a boundary con- 
dition consistent with hydrostatic pressure against the 
ice front, 


\4hvd u — ^-hA 
l 2 J x=L 




(17) 


while at the upstream end (x = 0) we set the velocity to 
be consistent with the imposed influx of ice. Advection 
of ice is determined from the continuity equation 


4. Coupled ice-ocean experiments 

Flaving seen that uncertainty in the poorly known 
drag coefficient can lead to a broad range of calculated 
melt rates, we will now examine the potential impact of 
this uncertainty on the flow of grounded ice. This 
investigation will require coupling the plume model 
to a reduced-dimensional model of ice stream and ice 
shelf flow. 

a. Ice model and coupling approach 

We use a version of the Dupont and Alley (2005) 
flowline model, which is a one-dimensional “shelfy stream” 
model (MacAyeal 1989). This type of model is the sim- 
plest that includes the longitudinal stresses responsible 
for ice shelf buttressing of inland ice flow. 

The momentum equation is given by 

d x (4hvd x u - = Pjghd x z b + h T y (u) + T b (u ) , 

(13) 

where h, 11 , and p, are the ice thickness, velocity, and 
density, respectively; Zb is the elevation of the ice base; 
and L y is the half width of the ice stream. The effective 
viscosity is defined as 

v^d x u\- y \ (14) 

where B = A~ 1/3 is the ice hardness parameter, and A is 
the ice softness parameter in Glen’s law. Lateral drag is 
parameterized by 

Ty = T y U V \ (15) 


d t h = —d x (uh) — m, (18) 

in which we neglect surface accumulation so that the 
basal melt rate m is the only forcing. We apply a constant 
flux boundary condition at the upstream end, while al- 
lowing free outflux of ice at the downstream boundary. 

The model is discretized using linear finite elements, 
with a Petrov-Galerkin upwinding method applied to 
(18). Nodal spacing is 250 m, with refinement to 10 m for 
at least 2 km up- and downstream of the grounding line 
when solving (13). This relatively high resolution is 
necessary for smooth grounding-line migration without 
numerical artifacts similar to those observed by Vieli 
and Payne (2005), despite our use of grounding-line in- 
terpolation and partial-element basal drag as in Parizek 
et al. (2010). Time discretization is fully implicit, with 24 
time steps per year. Coupling between ocean and ice 
models is handled as in Parizek and Walker (2010), with 
the ice shelf depth profile passed to the plume model and 
melt rates at ice shelf nodes returned. The full output of 
the plume model is also saved for analysis. 

Because the plume model is steady state, it can be 
called on time scales consistent with ice shelf evolution, 
avoiding the difficulty of allowing continuously changing 
domain geometry and the expense of running a time- 
dependent model at relatively fast ocean time scales. In 
contrast with earlier simplified coupled models (Walker 
and Flolland 2007; Walker et al. 2009), the time required 
to solve an ocean model consisting only of ordinary 
differential equations (ODEs) is negligible, allowing 
more computational resources to be dedicated to im- 
proved resolution of grounding-line migration. A simi- 
larly efficient approach is taken by Gladstone et al. 
(2012), who couple a flowline ice model with a box model 
of ocean circulation that is also a set of ODEs. In the 
experiments described below, the time between ocean 
model calls is one year. 


and basal drag by 


2206 


JOURNAL OF PHYSICAL OCEANOGRAPHY 


Volume 43 



0 50 100 150 


Along-flow distance (km) 

Fig. 4. Initial steady-state configuration for coupled model ex- 
periments shown by solid line. Final steady state of experiment 
with C d = 0.0035 and ocean temperature of — 0.6°C shown by 
dashed-dotted line. Final nonsteady configuration of experiment 
with C d = 0.0097 and ocean temperature of — 0.2°C just prior to 
grounding-line retreat out of model domain shown by dashed line. 

b. Assessing sensitivity to (ocean) drag coefficient 

In our coupled experiments, the plume model uses 
laboratory three-equation thermodynamics for each of 
the five values of the drag coefficient used in section 3; as 
discussed there, the uncertainty in melt rates due to this 
parameter can be much greater than that due to differ- 
ing thermodynamical schemes. To clearly assess the ef- 
fects of ocean parameters on the coupled system, it is 
necessary to begin with an ice shelf-ice stream that is 
a steady state of the ice model without oceanic forcing. 
Attempting to use the arbitrary ice shelf profiles given 
by (12) would introduce potentially large ice flow tran- 
sients, as seen by Parizek and Walker (2010), signifi- 
cantly complicating analysis. The initial ice configuration 
is identical to that used in Walker et al. (2008), with a 
roughly 99-km-long ice stream flowing up an inland- 
deepening bed (slope 3 X 10~ 3 ) and an ice shelf occu- 
pying the remainder of the 150-km domain (Fig. 4). The 
original study found that applying basal melting aver- 
aging 1 0 m yr 1 would cause grounding-line retreat to 
a new stable position if the melt was relatively evenly 
distributed across the shelf, but would cause unstable 
retreat (i.e., flotation of all ice in the domain) if the melt 
was sufficiently concentrated toward the grounding line. 
Basal melting averaging 15 m yr -1 was strong enough to 
cause unstable retreat in all cases. 

As might be expected, experiments with our coupled 
model at an ocean temperature of — 1.8°C show only 
minor differences in final steady state as C d is varied. 
Steady mean and maximum melt rates are 1.50 and 
1.79 myr -1 for C d = 0.0015, increasing to 2.43 and 



Fig. 5. Grounding-line retreat for coupled model with ocean 
temperature of — 0.6°C. 


2.94 myr -1 for C d = 0.0097. While this is a difference of 
over 60%, melt rates in all experiments remain small 
enough that grounding-line retreat ranges only from 
1.41 to 2.65 km. 

The effect of C d becomes slightly more pronounced as 
the ocean temperature is increased to — 1.2°C. Steady 
mean and maximum melt rates now range from 4.51 and 
4.74 myr -1 for C d = 0.0015 to 7.62 and 8.37 myr -1 for 
C d = 0.0097. Melt rates are still insufficient to cause 
complete flotation, but grounding-line retreat increases, 
ranging from 5.44 to 12.06 km. 

The most dramatic effect of varying C d is seen at an 
ocean temperature of — 0.6°C, where the drag coefficient 
determines stability (Fig. 5). For experiments reaching 
a new steady state, final mean and maximum melt rates 
range from 9.04 and 9.46 myr -1 at C d = 0.0015 to 11.79 
and 13.16 myr -1 at C d = 0.0035. While the latter run 
does have melt exceeding the 10 m yr~ 1 threshold at which 
unstable retreat was possible in Walker et al. (2008), its 
spatial distribution is only weakly concentrated near the 
grounding line, resulting in a long (41 km in 1260 years) but 
stable retreat. Both of the larger drag coefficient values 
introduced by Jenkins et al. (2010b) produce unstable re- 
treat, with C d = 0.0062 leading to flotation of the entire ice 
stream in 349 years, and C d = 0.0097 in 214 years. Mean 
and maximum melt rates at the end of these experiments 
are 14.18 and 17.77 myr -1 for C d = 0.0062, and 16.50 and 
20.80 myr -1 for C d = 0.0097. 

Finally, warming the ocean to — 0.2°C is sufficient to 
cause unstable retreat for all our values of C d . Final 
mean and maximum melt rates range from 13.52 and 
15.58 myr -1 for C d = 0.0015 up to 26.13 and 31.38myr _1 
for C d = 0.0097. This doubling of melt rate produces 
significant differences in the rate of retreat, with complete 
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Fig. 6. Max basal melt rates for coupled model with ocean 
temperature of — 0.6°C. 

flotation requiring as long as 554 years or as little as 95 
years. We note that the apparently slight difference be- 
tween C,i = 0.0015 and 0.0025 is enough to cause retreat 
times to vary by more than a factor of 2. 

In all experiments, melting increases rapidly in the 
first few decades as the ice shelf base steepens, demon- 
strating the feedback seen by Walker and Holland 
(2007). Once the ice stream has time to respond dy- 
namically to the onset of forcing, advection of ice mod- 
erates slopes and leads to melt rates increasing more 
slowly (for unstable runs) or reaching a steady state 
(Fig. 6). As the grounding line retreats, small variations 
(on the order of several centimeters per year) are su- 
perimposed on the overall trend in melt rate. When the 
ice at a model node thins to flotation, the slope of the 
first few kilometers of the ice shelf base increases slightly, 
leading to greater entrainment and a small but sudden 
increase in melting. As the newly floating ice adjusts to 
oceanic forcing, this increase fades over decades, until 
the ungrounding of another node restarts the cycle or 
a final grounding-line position is reached. Because the 
variations are generally less than 1% of the final melt 
rate, we expect that our asynchronous coupling method 
has not significantly affected our results. Rather, we point 
out subtle variations in calculated melt rates as evidence 
of the sensitivity of plume models to near-grounding-line 
slope. 

5. Conclusions 

The broad range of reduced-model experiments con- 
ducted in this study leads to three principal conclusions 
regarding ice-ocean thermodynamics. First, the con- 
stant Stanton number method proposed by Jenkins et al. 


(2010b) and the laboratory-based transfer velocities of 
Holland and Jenkins (1999) produce similar melt rates 
across a broad range of shelf geometries and ocean 
temperatures when both schemes are appropriately 
parameterized. The two-equation method (McPhee 1992) 
is also in reasonable agreement when basal slopes and 
ocean temperatures remain relatively low. The uncer- 
tainty in melt rates resulting from the choice of ther- 
modynamics can often be smaller than that resulting 
from uncertainty in a single parameter, the drag co- 
efficient at the ice-ocean interface, which influences the 
entrainment of warmer ambient water into the plume. 
Second, grounding-line retreat in a coupled ocean-ice 
shelf-ice stream model is highly sensitive to the drag 
coefficient, again because of the influence of this param- 
eter on entrainment. Third, as will be shown in the ap- 
pendix, ice shelves with very steep near-grounding-line 
slopes can experience lower melt rates than flatter shelves 
when high entrainment in this region limits buoyancy and 
slows the plume. 

Our use of reduced models allows exploration of a 
broad parameter space, but requires some caution when 
interpreting the results. While the plume model pro- 
duces a reasonable near-shelf boundary layer that cou- 
ples well with the interface thermodynamics, it also 
simplifies sub-ice shelf circulation by ignoring horizon- 
tal variations transverse to ice shelf flow and deep-ocean 
dynamics. The parameterization of entrainment, while 
based on observations, is simpler than the vertical mix- 
ing schemes of most ocean general circulation models, 
so that effects attributed here to the drag coefficient may 
be sensitive to multiple parameters in a more complex 
model. In addition, ice shelf grounding-line dynamics 
can be more complicated in nature than can be fully 
captured by a reduced-dimensional model on an ideal- 
ized domain. Nevertheless, our results strongly suggest 
that uncertainty about the effects of vertical mixing on 
ice-ocean thermodynamics limits the ability of coupled 
models to make accurate projections of grounding-line 
retreat. 

Given the scarcity of available data on sub-ice shelf 
ocean circulation, nearly any new observations would 
be useful for improving projections. Autonomous sub- 
mersibles (e.g., Jenkins et al. 2010a) can map bathyme- 
try and water properties, while the combination of 
phase-sensitive radar (Corr 2002; Jenkins et al. 2006) 
and thermistors deployed through boreholes (e.g., Nicholls 
et al. 2009) allows limited validation of parameteriza- 
tions of ice-ocean thermodynamics. However, as noted 
by Jenkins et al. (2010b), such observations constrain 
only the Stanton numbers and cannot be used to deter- 
mine the drag or turbulent exchange coefficients. Direct 
measurements of turbulent transfer at the ice shelf base 
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will be necessary to produce a clear picture of thermo- 
dynamical processes in the ice-ocean boundary layer 
and fully validate model parameterizations. Such mea- 
surements and their analysis are currently in progress for 
Pine Island, Larsen C (Nicholls et al. 2012), and George 
VI ice shelves, and the results, when available, are ex- 
pected to inform projections of sub-ice shelf ocean cir- 
culation and its role in forcing grounding-line retreat 
and sea level rise. 
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APPENDIX 

Effect of Steep Basal Slopes near the Grounding Line 

Although issues with plume separation prevent us 
from using the Holland et al. (2008) ice shelf profiles 
across a broad range of ocean temperature, we run some 
experiments with a deeper (1000 m) grounding line to 
examine the effect of relatively steep basal slopes. While 
high slopes near the grounding line are generally ex- 
pected to enhance melting, we find cases in which melt 
rates for the steepest ( n = 4) shelf are lower than melt 
rates for the n = 3 and 2 profiles. This results from near- 
grounding-line entrainment as well as the overall shape 
of the profiles, all of which have the same mean slope. 

These experiments span the full range of fixed ice 
shelf profiles described in section 3, with laboratory 
thermodynamics using all five previously considered 
values of C',/. For the highly concave (n = 4) profile, 
plume separation occurs within the first 50 km at lower 
ocean temperatures (<0°C for 275-km length and 
<-1.2°C for 550-km length) regardless of C,/. At higher 
ocean temperatures, at least one of the lower-exponent 
275-km-long shelves experiences greater melt than the 
highest-concavity shelf in all experiments. The 550-km 
shelves, which have lower concavity, display a slightly 
more complex dependence on ocean temperature and 
drag coefficient. For these profiles, once the ocean is warm 
enough to drive sufficient melting to prevent plume 
separation, there is a temperature range for which at 



Fig. Al. Melt rates for 275-km ice shelf profiles with deeper 
(1000 m) grounding line, with Cj = 0.0025 and ocean temperature 
of 0.8°C. 


least one lower-concavity shelf experiences greater melt 
than the highest-concavity shelf. (This range depends on 
the drag coefficient, extending 0.4°C warmer for C,/ = 
0.0015 than 0.0097.) When the ocean is warmed further, 
melting becomes sufficiently strong relative to entrain- 
ment that melt rates increase with concavity. 

As an example, we analyze runs with 275-km-long 
shelves, Cj = 0.0025, and ocean temperature of 0.8°C 
(Fig. Al). For this geometry, the n = 4 shelf has a max- 
imum slope of nearly 25°, roughly 3 times greater than 
for n = 3 (8.3°) and 10 times greater than for n = 2 (2.5°). 
The entrainment rate (5), which depends on the sine 
of the basal slope, follows roughly the same proportion 
in the first several hundred meters of these three ex- 
periments, resulting in a much thicker n = 4 plume at 

1 km along the flow (6.4 m, versus 2.6 m and 0.9 m, re- 
spectively). This close to the origin of the plume, melting 
has had little chance to contribute to buoyancy, so the 
momentum equation (2) requires a rapidly thickening 
plume to slow down; the n = 4 plume’s velocity falls 
below that of the n = 3 plume by 200 m along the flow 
and below that of the n = 2 plume by 900 m along the 
flow. Because the n = 4 profile remains steeper than the 
n = 3 and 2 profiles for less than 2.5 and 5 km, re- 
spectively, the n = 4 plume never fully recovers from 
this early deceleration, reaching a peak velocity of only 
17.3 cms -1 versus 21.6 and 22.9cms~ 1 for the n = 3 and 

2 plumes. Although the n = 2 plume reaches the highest 
velocity, it does so farther downstream and thus shal- 
lower than the n = 3 plume; the greater contrast at depth 
between the local freezing point and ambient ocean 
temperature results in the n = 3 plume producing the 
highest maximum melt rate. The dependence of melt 
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rate on plume velocity [through the turbulent exchange 
velocities (9)] then leads to the second- and third- 
steepest ice shelf profiles producing higher melt rates 
than the steepest. 

Should this result hold in nature, rather than being an 
artifact of a particular model, it presents the possibility 
of a limiting mechanism for melting of a coupled ocean- 
ice shelf system. With an early coupled model. Walker 
and Holland (2007) found a positive feedback between 
basal melting and ice shelf basal slope. At the relatively 
low slopes used in their experiments, it appeared that 
slope-melt feedback would grow indefinitely, until hal- 
ted by reaching equilibrium with ice flux into the shelf. 
The results of the present study suggest that in some 
cases the feedback could be slowed by the ice shelf base 
reaching a sufficiently steep slope, thus complicating the 
process by which the coupled system reaches equilib- 
rium, and possibly reducing grounding-line retreat. How- 
ever, we caution that this preliminary result requires 
confirmation by observations or by a three-dimensional 
ocean model capable of resolving the sub-shelf bound- 
ary layer. 
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